Session 6: Homework 3

1 Youth Risk Behavior Surveillance

Every two years, the Centers for Disease Control and Prevention conduct the Youth Risk Behavior Surveillance System (YRBSS) survey, where it takes data from high schoolers (9th through 12th grade), to analyze health patterns. You will work with a selected group of variables from a random sample of observations during one of the years the YRBSS was conducted.

1.1 Load the data

This data is part of the openintro textbook and we can load and inspect it. There are observations on 13 different variables, some categorical and some numerical. The meaning of each variable can be found by bringing up the help file:

?yrbss

data(yrbss)
glimpse(yrbss)

## Rows: 13,583
## Columns: 13
## $ age                      <int> 14, 14, 15, 15, 15, 15, 15, 14, 15, 15, 15, 1…
## $ gender                   <chr> "female", "female", "female", "female", "fema…
## $ grade                    <chr> "9", "9", "9", "9", "9", "9", "9", "9", "9", …
## $ hispanic                 <chr> "not", "not", "hispanic", "not", "not", "not"…
## $ race                     <chr> "Black or African American", "Black or Africa…
## $ height                   <dbl> NA, NA, 1.73, 1.60, 1.50, 1.57, 1.65, 1.88, 1…
## $ weight                   <dbl> NA, NA, 84.4, 55.8, 46.7, 67.1, 131.5, 71.2, …
## $ helmet_12m               <chr> "never", "never", "never", "never", "did not …
## $ text_while_driving_30d   <chr> "0", NA, "30", "0", "did not drive", "did not…
## $ physically_active_7d     <int> 4, 2, 7, 0, 2, 1, 4, 4, 5, 0, 0, 0, 4, 7, 7, …
## $ hours_tv_per_school_day  <chr> "5+", "5+", "5+", "2", "3", "5+", "5+", "5+",…
## $ strength_training_7d     <int> 0, 0, 0, 0, 1, 0, 2, 0, 3, 0, 3, 0, 0, 7, 7, …
## $ school_night_hours_sleep <chr> "8", "6", "<5", "6", "9", "8", "9", "6", "<5"…

Before you carry on with your analysis, it’s is always a good idea to check with skimr::skim() to get a feel for missing values, summary statistics of numerical variables, and a very rough histogram.

1.2 Exploratory Data Analysis

You will first start with analyzing the weight of participants in kilograms. Using visualization and summary statistics, describe the distribution of weights. How many observations are we missing weights from?

# find missing values and visualize
  skim(yrbss)

Data summary

Name	yrbss
Number of rows	13583
Number of columns	13
_______________________
Column type frequency:
character	8
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
gender	12	1.00	4	6	0	2	0
grade	79	0.99	1	5	0	5	0
hispanic	231	0.98	3	8	0	2	0
race	2805	0.79	5	41	0	5	0
helmet_12m	311	0.98	5	12	0	6	0
text_while_driving_30d	918	0.93	1	13	0	8	0
hours_tv_per_school_day	338	0.98	1	12	0	7	0
school_night_hours_sleep	1248	0.91	1	3	0	7	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
age	77	0.99	16.16	1.26	12.00	15.0	16.00	17.00	18.00	▁▂▅▅▇
height	1004	0.93	1.69	0.10	1.27	1.6	1.68	1.78	2.11	▁▅▇▃▁
weight	1004	0.93	67.91	16.90	29.94	56.2	64.41	76.20	180.99	▆▇▂▁▁
physically_active_7d	273	0.98	3.90	2.56	0.00	2.0	4.00	7.00	7.00	▆▂▅▃▇
strength_training_7d	1176	0.91	2.95	2.58	0.00	0.0	3.00	5.00	7.00	▇▂▅▂▅

  yrbss %>% 
    ggplot(aes(x = weight))+
    geom_density()

Missing 1004 observations.

A more insightful graph would be something like this:

yrbss %>% 
  filter (!is.na(physically_active_7d)) %>% 
ggplot(aes(x = factor(physically_active_7d))) + 
  geom_boxplot(aes(y = weight)) +
  labs(x="Activity level",
       y="weight")

Next, consider the possible relationship between a high schooler’s weight and their physical activity. Plotting the data is a useful first step because it helps us quickly visualize trends, identify strong associations, and develop research questions.

Let’s create a new variable in the dataframe yrbss, called physical_3plus , which will be yes if they are physically active for at least 3 days a week, and no otherwise. You may also want to calculate the number and % of those who are and are not active for more than 3 days. Use the count() function and see if you get the same results as group_by()... summarise()

# create a new variable
  yrbss<- yrbss %>% 
  mutate(physical_3plus = if_else(physically_active_7d >= 3, "yes", "no")) %>% 
  drop_na(physical_3plus)

# number and percentage using count()
  physical_act_matrix <- yrbss %>% 
    count(physical_3plus) %>% 
    summarise(count = n, 
            percentage = n/sum(n))

# number and percentage using group_by()
  physical_act_matrix_v2 <- yrbss %>% 
   group_by(physical_3plus) %>% 
   summarise(count = n())  %>% 
   mutate(percentage = count/sum(count))

Can you provide a 95% confidence interval for the population proportion of high schools that are NOT active 3 or more days per week?

ci_prop_not_excercise <- physical_act_matrix_v2 %>% 
  group_by(physical_3plus) %>% 
  summarize(mean = percentage,
            se = sqrt(percentage * (1 - percentage) / count),
            margin = se * qnorm(0.975),
            upper = mean + margin,
            lower = mean - margin
            )

ci_prop_not_excercise

## # A tibble: 2 × 6
##   physical_3plus  mean      se  margin upper lower
##   <chr>          <dbl>   <dbl>   <dbl> <dbl> <dbl>
## 1 no             0.331 0.00709 0.0139  0.345 0.317
## 2 yes            0.669 0.00499 0.00977 0.679 0.659

#interval is (0.3003; 0.3357)

Make a boxplot of physical_3plus vs. weight. Is there a relationship between these two variables? What did you expect and why?

ggplot(yrbss, aes (x = weight, y = physical_3plus)) +
  geom_boxplot() +
  theme_bw() +
  labs(x = "weight",
       y = "Is person physically active?")

No obvious relationship can be spotted from this particular graph. Didn’t expect anything since we should’t bring our biases into the classroom.

1.3 Confidence Interval

Boxplots show how the medians of the two distributions compare, but we can also compare the means of the distributions using either a confidence interval or a hypothesis test. Note that when we calculate the mean, SD, etc. weight in these groups using the mean function, we must ignore any missing values by setting the na.rm = TRUE.

ci <- yrbss  %>% 
  drop_na() %>% 
  group_by(physical_3plus) %>% 
  summarise(mean = mean(weight),
            min = min(weight),
            max = max(weight),
            median = median(weight),
            stdev = sd(weight),
            se = stdev/sqrt(n()),
            t_crit = qt(0.975, n() - 1 ),
            upper = mean + t_crit * se,
            lower = mean - t_crit * se) 
 
ci %>% 
  select(physical_3plus, mean, lower, upper) 

## # A tibble: 2 × 4
##   physical_3plus  mean lower upper
##   <chr>          <dbl> <dbl> <dbl>
## 1 no              67.1  66.5  67.8
## 2 yes             68.7  68.3  69.1

CIs dont overlap sowe can infer statistically significant difference.

1.4 Hypothesis test with formula

Write the null and alternative hypotheses for testing whether mean weights are different for those who exercise at least times a week and those who don’t.

t.test(weight ~ physical_3plus, data = yrbss)

## 
##  Welch Two Sample t-test
## 
## data:  weight by physical_3plus
## t = -5, df = 7479, p-value = 9e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.42 -1.12
## sample estimates:
##  mean in group no mean in group yes 
##              66.7              68.4

1.5 Hypothesis test with infer

Next, we will introduce a new function, hypothesize, that falls into the infer workflow. You will use this method for conducting hypothesis tests.

But first, we need to initialize the test, which we will save as obs_diff.

obs_diff <- yrbss %>%
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

Notice how you can use the functions specify and calculate again like you did for calculating confidence intervals. Here, though, the statistic you are searching for is the difference in means, with the order being yes - no != 0.

After you have initialized the test, you need to simulate the test on the null distribution, which we will save as null.

null_dist <- yrbss %>%
  # specify variables
  specify(weight ~ physical_3plus) %>%
  
  # assume independence, i.e, there is no difference
  hypothesize(null = "independence") %>%
  
  # generate 1000 reps, of type "permute"
  generate(reps = 1000, type = "permute") %>%
  
  # calculate statistic of difference, namely "diff in means"
  calculate(stat = "diff in means", order = c("yes", "no"))

Here, hypothesize is used to set the null hypothesis as a test for independence, i.e., that there is no difference between the two population means. In one sample cases, the null argument can be set to point to test a hypothesis relative to a point estimate.

Also, note that the type argument within generate is set to permute, which is the argument when generating a null distribution for a hypothesis test.

We can visualize this null distribution with the following code:

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram()

Now that the test is initialized and the null distribution formed, we can visualise to see how many of these null permutations have a difference of at least obs_stat of 1.77?

We can also calculate the p-value for your hypothesis test using the function infer::get_p_value().

null_dist %>% visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "two-sided")

null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "two_sided")

## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

This the standard workflow for performing hypothesis tests.

2 IMDB ratings: Differences between directors

Recall the IMBD ratings data. I would like you to explore whether the mean IMDB rating for Steven Spielberg and Tim Burton are the same or not. I have already calculated the confidence intervals for the mean ratings of these two directors and as you can see they overlap.

First, I would like you to reproduce this graph. You may find geom_errorbar() and geom_rect() useful.

In addition, you will run a hpothesis test. You should use both the t.test command and the infer package to simulate from a null distribution, where you assume zero difference between the two.

Before anything, write down the null and alternative hypotheses, as well as the resulting test statistic and the associated t-stat or p-value. At the end of the day, what do you conclude?

You can load the data and examine its structure

movies <- read_csv(here::here("data", "movies.csv"))
glimpse(movies)

## Rows: 2,961
## Columns: 11
## $ title               <chr> "Avatar", "Titanic", "Jurassic World", "The Avenge…
## $ genre               <chr> "Action", "Drama", "Action", "Action", "Action", "…
## $ director            <chr> "James Cameron", "James Cameron", "Colin Trevorrow…
## $ year                <dbl> 2009, 1997, 2015, 2012, 2008, 1999, 1977, 2015, 20…
## $ duration            <dbl> 178, 194, 124, 173, 152, 136, 125, 141, 164, 93, 1…
## $ gross               <dbl> 7.61e+08, 6.59e+08, 6.52e+08, 6.23e+08, 5.33e+08, …
## $ budget              <dbl> 2.37e+08, 2.00e+08, 1.50e+08, 2.20e+08, 1.85e+08, …
## $ cast_facebook_likes <dbl> 4834, 45223, 8458, 87697, 57802, 37723, 13485, 920…
## $ votes               <dbl> 886204, 793059, 418214, 995415, 1676169, 534658, 9…
## $ reviews             <dbl> 3777, 2843, 1934, 2425, 5312, 3917, 1752, 1752, 35…
## $ rating              <dbl> 7.9, 7.7, 7.0, 8.1, 9.0, 6.5, 8.7, 7.5, 8.5, 7.2, …

Your R code and analysis should go here. If you want to insert a blank chunk of R code you can just hit Ctrl/Cmd+Alt+I

#obtaining confidence interval as a pre condition for replicating the graph

movies_Spielberg_burton <- movies %>%
  group_by(director) %>%
  filter(director %in% c("Steven Spielberg" , "Tim Burton"))

#glimpse(movies_Spielberg_burton)
formula_ci <- movies_Spielberg_burton %>%
summarize(mean_rating = mean(rating, na.rm = TRUE),
                      sd_rating = sd(rating, na.rm=TRUE),
                      count = n(),
                      se_rating = sd_rating / sqrt(count),
                      lower_cl = mean_rating - qt(0.975, count - 1) * se_rating,
                      upper_cl = mean_rating + qt(0.975, count - 1) * se_rating
                      )
formula_ci

## # A tibble: 2 × 7
##   director         mean_rating sd_rating count se_rating lower_cl upper_cl
##   <chr>                  <dbl>     <dbl> <int>     <dbl>    <dbl>    <dbl>
## 1 Steven Spielberg        7.57     0.695    23     0.145     7.27     7.87
## 2 Tim Burton              6.93     0.749    16     0.187     6.53     7.33

#graph
  ggplot(formula_ci, aes(x = mean_rating,
                             y = director, 
                             color = director)) +
          
  #create errorbar
  geom_errorbar(aes(xmin = lower_cl, xmax = upper_cl, color = director),
                width = 0.1,
                size =1.5,
                show.legend = FALSE) +
  
  #plot a dot to show the mean on the errorbar
  stat_summary(fun.y = mean, geom="point", size = 4, show.legend = FALSE) +
  
  # add label to the upper end of the errorbar
  geom_text(aes(label = round(upper_cl,2), x = upper_cl),
            vjust = -1,colour = "black") +
  
  
  # add label to the lower end of the errorbar
  geom_text(aes(label = round(lower_cl,2), x = lower_cl),
            vjust = -1,colour = "black") +
  
  #add label to the "mean dot" on the errorbar
  geom_text(aes(label = round(mean_rating,2)),
            vjust = -1,size =6,colour = "black") +
  
  #use geom_rect to show the overlap
  geom_rect(aes(xmin = min(upper_cl),
                xmax=max(lower_cl),
                ymin=-Inf,
                ymax=Inf),
            fill = "grey",
            color= "grey",
            alpha = 0.3) +
  theme_bw() +
  labs(title="Do Spielberg and Burton have the same IMDB rating",
       subtitle="95% confidence interval overlap",
       x="Mean IMDB rating",
       y="") +
  NULL

#Hypothesis test

#Ho=mean difference equals to 0 #H1=mean difference does not equal to 0

  t.test(rating~director, data=movies_Spielberg_burton)

## 
##  Welch Two Sample t-test
## 
## data:  rating by director
## t = 3, df = 31, p-value = 0.01
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.16 1.13
## sample estimates:
## mean in group Steven Spielberg       mean in group Tim Burton 
##                           7.57                           6.93

  #library(infer)
set.seed(1234)
ratings_spielberg_burton <- movies_Spielberg_burton %>%
  specify(rating~director) %>%
  hypothesize(null="independence") %>%
  generate(reps=1000, type="permute") %>%
  calculate(stat = "diff in means", order=c("Steven Spielberg","Tim Burton"))

  ratings_spielberg_burton %>% visualize()

3 Omega Group plc- Pay Discrimination

At the last board meeting of Omega Group Plc., the headquarters of a large multinational company, the issue was raised that women were being discriminated in the company, in the sense that the salaries were not the same for male and female executives. A quick analysis of a sample of 50 employees (of which 24 men and 26 women) revealed that the average salary for men was about 8,700 higher than for women. This seemed like a considerable difference, so it was decided that a further analysis of the company salaries was warranted.

You are asked to carry out the analysis. The objective is to find out whether there is indeed a significant difference between the salaries of men and women, and whether the difference is due to discrimination or whether it is based on another, possibly valid, determining factor.

3.1 Loading the data

omega <- read_csv(here::here("data", "omega.csv"))
glimpse(omega) # examine the data frame

## Rows: 50
## Columns: 3
## $ salary     <dbl> 81894, 69517, 68589, 74881, 65598, 76840, 78800, 70033, 635…
## $ gender     <chr> "male", "male", "male", "male", "male", "male", "male", "ma…
## $ experience <dbl> 16, 25, 15, 33, 16, 19, 32, 34, 1, 44, 7, 14, 33, 19, 24, 3…

3.2 Relationship Salary - Gender ?

The data frame omega contains the salaries for the sample of 50 executives in the company. Can you conclude that there is a significant difference between the salaries of the male and female executives?

Note that you can perform different types of analyses, and check whether they all lead to the same conclusion

. Confidence intervals . Hypothesis testing . Correlation analysis . Regression

Calculate summary statistics on salary by gender. Also, create and print a dataframe where, for each gender, you show the mean, SD, sample size, the t-critical, the SE, the margin of error, and the low/high endpoints of a 95% condifence interval

# Summary Statistics of salary by gender
mosaic::favstats (salary ~ gender, data=omega)

##   gender   min    Q1 median    Q3   max  mean   sd  n missing
## 1 female 47033 60338  64618 70033 78800 64543 7567 26       0
## 2   male 54768 68331  74675 78568 84576 73239 7463 24       0

# Summary Statistics of salary by gender
df<-mosaic::favstats (salary ~ gender, data=omega)
df

##   gender   min    Q1 median    Q3   max  mean   sd  n missing
## 1 female 47033 60338  64618 70033 78800 64543 7567 26       0
## 2   male 54768 68331  74675 78568 84576 73239 7463 24       0

# Dataframe with two rows (male-female) and having as columns gender, mean, SD, sample size, 

df_01<-select(df,gender,mean,sd,n)

# the t-critical value, the standard error, the margin of error, 

df_02 <- df_01 %>% 
  mutate(t_critical = qt(0.975, n-1),
         SE = sd / sqrt(n),
         margin_error = SE * t_critical)
df_02

##   gender  mean   sd  n t_critical   SE margin_error
## 1 female 64543 7567 26       2.06 1484         3056
## 2   male 73239 7463 24       2.07 1523         3151

# and the low/high endpoints of a 95% condifence interval

df_03 <- df_02 %>% 
  mutate(high_end = mean + margin_error,
         low_end = mean - margin_error)
df_03

##   gender  mean   sd  n t_critical   SE margin_error high_end low_end
## 1 female 64543 7567 26       2.06 1484         3056    67599   61486
## 2   male 73239 7463 24       2.07 1523         3151    76390   70088

What can you conclude from your analysis? A couple of sentences would be enough

There seems to be a significant difference between female and male salaries. In fact, while the confidence interval for females’ salaries goes from 61,486 to 67,599, the one for men goes from 70,088 to 76,390. In other words, the two CIs do not overlap, something that suggests us to directly reject the null hypothesis that the two means are equal.

You can also run a hypothesis testing, assuming as a null hypothesis that the mean difference in salaries is zero, or that, on average, men and women make the same amount of money. You should tun your hypothesis testing using t.test() and with the simulation method from the infer package.

# hypothesis testing using t.test() 

t.test(salary~gender,data=omega)

## 
##  Welch Two Sample t-test
## 
## data:  salary by gender
## t = -4, df = 48, p-value = 2e-04
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -12973  -4420
## sample estimates:
## mean in group female   mean in group male 
##                64543                73239

glimpse(omega)

## Rows: 50
## Columns: 3
## $ salary     <dbl> 81894, 69517, 68589, 74881, 65598, 76840, 78800, 70033, 635…
## $ gender     <chr> "male", "male", "male", "male", "male", "male", "male", "ma…
## $ experience <dbl> 16, 25, 15, 33, 16, 19, 32, 34, 1, 44, 7, 14, 33, 19, 24, 3…

# hypothesis testing using infer package

 diff_salary <- omega %>%
   specify(salary ~ gender) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat ="diff in means",
             order=c("male","female"))
 diff_salary %>% 
   visualize()

What can you conclude from your analysis? A couple of sentences would be enough

First of all, since the p value is lower than 0.05, we reject the null hypothesis at a 95% confidence level. In other words, we can say that there is a statistically significant difference between the means man and female salaries.

3.3 Relationship Experience - Gender?

At the board meeting, someone raised the issue that there was indeed a substantial difference between male and female salaries, but that this was attributable to other reasons such as differences in experience. A questionnaire send out to the 50 executives in the sample reveals that the average experience of the men is approximately 21 years, whereas the women only have about 7 years experience on average (see table below).

# Summary Statistics of experience by gender
favstats (experience ~ gender, data=omega)

##   gender min    Q1 median   Q3 max  mean    sd  n missing
## 1 female   0  0.25    3.0 14.0  29  7.38  8.51 26       0
## 2   male   1 15.75   19.5 31.2  44 21.12 10.92 24       0

Based on this evidence, can you conclude that there is a significant difference between the experience of the male and female executives? Perform similar analyses as in the previous section. Does your conclusion validate or endanger your conclusion about the difference in male and female salaries?

# Summary Statistics of experience by gender
df<-mosaic::favstats (experience ~ gender, data=omega)
df

##   gender min    Q1 median   Q3 max  mean    sd  n missing
## 1 female   0  0.25    3.0 14.0  29  7.38  8.51 26       0
## 2   male   1 15.75   19.5 31.2  44 21.12 10.92 24       0

# Dataframe with two rows (male-female) and having as columns gender, mean, SD, sample size, 

df_05<-select(df,gender,mean,sd,n)

# the t-critical value, the standard error, the margin of error, 

df_06<-df_05 %>% 
  mutate(t_critical=qt(0.975,n-1),
         SE=sd/sqrt(n),
         margin_error=SE*t_critical)
df_06

##   gender  mean    sd  n t_critical   SE margin_error
## 1 female  7.38  8.51 26       2.06 1.67         3.44
## 2   male 21.12 10.92 24       2.07 2.23         4.61

# and the low/high endpoints of a 95% condifence interval

df_07<-df_06 %>% 
  mutate(high_end=mean+margin_error,
         low_end=mean-margin_error)
df_07

##   gender  mean    sd  n t_critical   SE margin_error high_end low_end
## 1 female  7.38  8.51 26       2.06 1.67         3.44     10.8    3.95
## 2   male 21.12 10.92 24       2.07 2.23         4.61     25.7   16.52

# hypothesis testing using t.test() 

t.test(experience~gender,data=omega)

## 
##  Welch Two Sample t-test
## 
## data:  experience by gender
## t = -5, df = 43, p-value = 1e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -19.35  -8.13
## sample estimates:
## mean in group female   mean in group male 
##                 7.38                21.12

glimpse(omega)

## Rows: 50
## Columns: 3
## $ salary     <dbl> 81894, 69517, 68589, 74881, 65598, 76840, 78800, 70033, 635…
## $ gender     <chr> "male", "male", "male", "male", "male", "male", "male", "ma…
## $ experience <dbl> 16, 25, 15, 33, 16, 19, 32, 34, 1, 44, 7, 14, 33, 19, 24, 3…

# hypothesis testing using infer package

 diff_experience <- omega %>%
   specify(experience ~ gender) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat ="diff in means",
             order=c("male","female"))
 diff_experience %>% 
   visualize()

At any reasonable significance level we can reject the H0 in favour of H1 - there is a significant difference between the experience of executives of different gender. This can be an explanation for salary difference.

3.4 Relationship Salary - Experience ?

Someone at the meeting argues that clearly, a more thorough analysis of the relationship between salary and experience is required before any conclusion can be drawn about whether there is any gender-based salary discrimination in the company.

Analyse the relationship between salary and experience. Draw a scatterplot to visually inspect the data

omega %>% 
ggplot(aes(x = experience, y = salary, color = gender))+
  geom_point() +
  geom_smooth(alpha=0.3) +
  theme_bw() +
  labs(x = "Experience, Years",
       y = "Salary",
       title = "Relationship between experience and salary") +
  theme(legend.position = "none") +
  facet_wrap(~ gender, scales = "free")

3.5 Check correlations between the data

You can use GGally:ggpairs() to create a scatterplot and correlation matrix. Essentially, we change the order our variables will appear in and have the dependent variable (Y), salary, as last in our list. We then pipe the dataframe to ggpairs() with aes arguments to colour by gender and make ths plots somewhat transparent (alpha = 0.3).

omega %>% 
  select(gender, experience, salary) %>% #order variables they will appear in ggpairs()
  ggpairs(aes (colour = gender, alpha = 0.3))+
  theme_bw()

Look at the salary vs experience scatterplot. What can you infer from this plot? Explain in a couple of sentences

From the plot we can infer that salary and experience are positively correlated. In fact, salary increases as experience increases. We can also infer that salary variability is highest in the first years of experiences, then it decreases when experiences in between 1 and 15 years, and then it increases again for employees with several years of experience.

4 Challenge 1: Yield Curve inversion

Every so often, we hear warnings from commentators on the “inverted yield curve” and its predictive power with respect to recessions. An explainer what a inverted yield curve is can be found here. If you’d rather listen to something, here is a great podcast from NPR on yield curve indicators

In addition, many articles and commentators think that, e.g., Yield curve inversion is viewed as a harbinger of recession. One can always doubt whether inversions are truly a harbinger of recessions, and use the attached parable on yield curve inversions.

In our case we will look at US data and use the FRED database to download historical yield curve rates, and plot the yield curves since 1999 to see when the yield curves flatten. If you want to know more, a very nice article that explains the yield curve is and its inversion can be found here. At the end of this challenge you should produce this chart

First, we will load the yield curve data file that contains data on the yield curve since 1960-01-01

yield_curve <- read_csv(here::here("data", "yield_curve.csv"))

glimpse(yield_curve)

## Rows: 6,884
## Columns: 5
## $ date      <date> 1960-01-01, 1960-02-01, 1960-03-01, 1960-04-01, 1960-05-01,…
## $ series_id <chr> "TB3MS", "TB3MS", "TB3MS", "TB3MS", "TB3MS", "TB3MS", "TB3MS…
## $ value     <dbl> 4.35, 3.96, 3.31, 3.23, 3.29, 2.46, 2.30, 2.30, 2.48, 2.30, …
## $ maturity  <chr> "3m", "3m", "3m", "3m", "3m", "3m", "3m", "3m", "3m", "3m", …
## $ duration  <chr> "3-Month Treasury Bill", "3-Month Treasury Bill", "3-Month T…

Our dataframe yield_curve has five columns (variables):

• date: already a date object
• series_id: the FRED database ticker symbol
• value: the actual yield on that date
• maturity: a short hand for the maturity of the bond
• duration: the duration, written out in all its glory!

4.1 Plotting the yield curve

This may seem long but it should be easy to produce the following three plots

4.1.1 Yields on US rates by duration since 1960

4.1.2 Monthly yields on US rates by duration since 1999 on a year-by-year basis

4.1.3 3-month and 10-year yields since 1999

According to Wikipedia’s list of recession in the United States, since 1999 there have been two recession in the US: between Mar 2001–Nov 2001 and between Dec 2007–June 2009. Does the yield curve seem to flatten before these recessions? Can a yield curve flattening really mean a recession is coming in the US? Since 1999, when did short-term (3 months) yield more than longer term (10 years) debt?

Besides calculating the spread (10year - 3months), there are a few things we need to do to produce our final plot

1. Setup data for US recessions
2. Superimpose recessions as the grey areas in our plot
3. Plot the spread between 30 years and 3 months as a blue/red ribbon, based on whether the spread is positive (blue) or negative(red)

• For the first, the code below creates a dataframe with all US recessions since 1946

# get US recession dates after 1946 from Wikipedia 
# https://en.wikipedia.org/wiki/List_of_recessions_in_the_United_States

recessions <- tibble(
  from = c("1948-11-01", "1953-07-01", "1957-08-01", "1960-04-01", "1969-12-01", "1973-11-01", "1980-01-01","1981-07-01", "1990-07-01", "2001-03-01", "2007-12-01","2020-02-01"),  
  to = c("1949-10-01", "1954-05-01", "1958-04-01", "1961-02-01", "1970-11-01", "1975-03-01", "1980-07-01", "1982-11-01", "1991-03-01", "2001-11-01", "2009-06-01", "2020-04-30") 
  )  %>% 
  mutate(From = ymd(from), 
         To=ymd(to),
         duration_days = To-From)


recessions

## # A tibble: 12 × 5
##    from       to         From       To         duration_days
##    <chr>      <chr>      <date>     <date>     <drtn>       
##  1 1948-11-01 1949-10-01 1948-11-01 1949-10-01 334 days     
##  2 1953-07-01 1954-05-01 1953-07-01 1954-05-01 304 days     
##  3 1957-08-01 1958-04-01 1957-08-01 1958-04-01 243 days     
##  4 1960-04-01 1961-02-01 1960-04-01 1961-02-01 306 days     
##  5 1969-12-01 1970-11-01 1969-12-01 1970-11-01 335 days     
##  6 1973-11-01 1975-03-01 1973-11-01 1975-03-01 485 days     
##  7 1980-01-01 1980-07-01 1980-01-01 1980-07-01 182 days     
##  8 1981-07-01 1982-11-01 1981-07-01 1982-11-01 488 days     
##  9 1990-07-01 1991-03-01 1990-07-01 1991-03-01 243 days     
## 10 2001-03-01 2001-11-01 2001-03-01 2001-11-01 245 days     
## 11 2007-12-01 2009-06-01 2007-12-01 2009-06-01 548 days     
## 12 2020-02-01 2020-04-30 2020-02-01 2020-04-30  89 days

recessions_ <- recessions %>% 
  filter(year(From) >=1959)

#calculate yield differences
yield_curve_graph <- yield_curve %>% 
  select(date, duration, value) %>% 
  pivot_wider(names_from = duration, values_from = value) %>% 
  janitor::clean_names() %>% 
  mutate(diff = x10_year_treasury_rate - x3_month_treasury_bill)

#add recessions
data_for_last_challenge_1_graph_xd <- merge(yield_curve_graph, recessions_)


ggplot(data_for_last_challenge_1_graph_xd, aes(date, diff)) +
  geom_line() +
  geom_line(aes(y = 0), 
            color = "black") +
  geom_rect(aes(xmin = From, 
                xmax = To, 
                ymin = -3, 
                ymax=5), 
            fill="grey") +
  geom_ribbon(aes(ymin = pmin(diff, 0), 
                  ymax = 0), 
              fill = "red", 
              alpha = 0.3) +
  geom_ribbon(aes(ymin = 0, 
                  ymax = pmax(diff, 0)), 
              fill = "blue", 
              alpha = 0.3) +
  geom_rug(sides = 'b', 
           data = subset(yield_curve_graph,
                         diff >= 0 ), 
           color = "blue", 
           alpha = 0.3) +
  geom_rug(sides = 'b', 
           data= subset(yield_curve_graph, 
                        diff < 0 ), 
           color = "red", 
           alpha = 0.3) +
  theme(legend.position = "none") +
  theme_minimal() +
  labs (title = "Yield Curve Inversion: 10-year minus 3-month U.S.Treasury rates",
        subtitle = "difference in %, shaded area corresponds to recession",
        caption = "Sources: St Louis Federal Reserve Economic Database (FRED)",
        y = "Difference (10 year - 3 month) yield in %",
        x = "") 

• To add the grey shaded areas corresponding to recessions, we use geom_rect()
• to colour the ribbons blue/red we must see whether the spread is positive or negative and then use geom_ribbon(). You should be familiar with this from last week’s homework on the excess weekly/monthly rentals of Santander Bikes in London.

5 Challenge 2: GDP components over time and among countries

At the risk of oversimplifying things, the main components of gross domestic product, GDP are personal consumption (C), business investment (I), government spending (G) and net exports (exports - imports). You can read more about GDP and the different approaches in calculating at the Wikipedia GDP page.

The GDP data we will look at is from the United Nations’ National Accounts Main Aggregates Database, which contains estimates of total GDP and its components for all countries from 1970 to today. We will look at how GDP and its components have changed over time, and compare different countries and how much each component contributes to that country’s GDP. The file we will work with is GDP and its breakdown at constant 2010 prices in US Dollars and it has already been saved in the Data directory. Have a look at the Excel file to see how it is structured and organised

UN_GDP_data  <-  read_excel(here::here('data',"Download-GDPconstant-USD-countries.xls"), # Excel filename, alter it to make sure it match path on indicidual's computer
                sheet="Download-GDPconstant-USD-countr", # Sheet name
                skip=2) # Number of rows to skip

The first thing you need to do is to tidy the data, as it is in wide format and you must make it into long, tidy format. Please express all figures in billions (divide values by 1e9, or \(10^9\)), and you want to rename the indicators into something shorter.

tidy_GDP_data  <-  UN_GDP_data %>% 
   pivot_longer(cols=c(4:51), 
               names_to = "year", 
               values_to = "gdp", 
               values_drop_na = FALSE) %>% 
  mutate(gdp_billion=gdp/1e9) %>% #express all figures in billions
  select(c(1:4,6)) %>% 
  pivot_wider(names_from=IndicatorName,
              values_from=gdp_billion) %>% # to rename the indicator
  
  rename("Exports"="Exports of goods and services",
         "Imports"="Imports of goods and services",
         "Government expenditure"="General government final consumption expenditure",
         "Household expenditure"="Household consumption expenditure (including Non-profit institutions serving households)") 
  
  
glimpse(tidy_GDP_data)

## Rows: 10,560
## Columns: 20
## $ CountryID                                                                            <dbl> …
## $ Country                                                                              <chr> …
## $ year                                                                                 <chr> …
## $ `Final consumption expenditure`                                                      <dbl> …
## $ `Household expenditure`                                                              <dbl> …
## $ `Government expenditure`                                                             <dbl> …
## $ `Gross capital formation`                                                            <dbl> …
## $ `Gross fixed capital formation (including Acquisitions less disposals of valuables)` <dbl> …
## $ Exports                                                                              <dbl> …
## $ Imports                                                                              <dbl> …
## $ `Gross Domestic Product (GDP)`                                                       <dbl> …
## $ `Agriculture, hunting, forestry, fishing (ISIC A-B)`                                 <dbl> …
## $ `Mining, Manufacturing, Utilities (ISIC C-E)`                                        <dbl> …
## $ `Manufacturing (ISIC D)`                                                             <dbl> …
## $ `Construction (ISIC F)`                                                              <dbl> …
## $ `Wholesale, retail trade, restaurants and hotels (ISIC G-H)`                         <dbl> …
## $ `Transport, storage and communication (ISIC I)`                                      <dbl> …
## $ `Other Activities (ISIC J-P)`                                                        <dbl> …
## $ `Total Value Added`                                                                  <dbl> …
## $ `Changes in inventories`                                                             <dbl> …

# Let us compare GDP components for these 3 countries
country_list <- c("United States","India", "Germany")

selected_data<-tidy_GDP_data %>% 
  #select the data we need
  filter(Country %in% country_list) %>% 
  select(c("Country","year",
           "Gross capital formation",
           "Exports","Imports",
           "Government expenditure",
           "Household expenditure")) %>% 
  
  #in order to color the data, we have to pivot_longer first
  pivot_longer(cols=c("Gross capital formation","Exports","Imports",
           "Government expenditure","Household expenditure"),
              names_to="Components of GDP",values_to="gdp_per_indicator") 

selected_data %>% 
  group_by(`Components of GDP`) %>% 
  ggplot(aes(x=as.numeric(year),y=gdp_per_indicator,color=`Components of GDP`))+
  geom_line(size=1.1)+
  theme_bw()+
  labs(title="GDP components over time",
       subtitle="In constant 2010 USD",
       x=" ",
       y="Billion US$")+
  scale_x_continuous(breaks=seq(1970,2020,by=10))+
  facet_wrap(~Country)

First, can you produce this plot?

GDP_proportion <- tidy_GDP_data %>% 
  mutate(GDP = `Gross capital formation` + `Exports` - `Imports` +
           `Government expenditure` + `Household expenditure`,
         "Household Expenditure" = `Household expenditure`/GDP,
         "Gross Capital Formation" = `Gross capital formation`/GDP,
         "Government Expenditure" = `Government expenditure`/GDP,
         "Net Exports" = (`Exports`-`Imports`)/GDP) %>% 
  select(c(1:3,22:25)) %>% 
  filter(Country %in% c("China", "India", "United States")) %>% 
  pivot_longer(cols = c("Household Expenditure", "Gross Capital Formation",
                      "Government Expenditure", "Net Exports"),
              names_to = "Components of GDP", values_to = "proportion") 

GDP_proportion %>% 
  group_by(`Components of GDP`) %>% 
  ggplot(aes(x = as.numeric(year), y = proportion,color = `Components of GDP`))+
  geom_line(size = 1.1)+
  theme_bw()+
  labs(title = "GDP and its breakdown at constant2010 prices in US Dollars",
       caption = "Source: United Nations,
       https://unstates.un.org/unsd/snaama/Downloads",
       x = " ",
       color = " ")+
  scale_x_continuous(breaks = seq(1970, 2020, by = 10))+
  facet_wrap(~Country)

Secondly, recall that GDP is the sum of Household Expenditure (Consumption C), Gross Capital Formation (business investment I), Government Expenditure (G) and Net Exports (exports - imports). Even though there is an indicator Gross Domestic Product (GDP) in your dataframe, I would like you to calculate it given its components discussed above.

What is the % difference between what you calculated as GDP and the GDP figure included in the dataframe?

What is this last chart telling you? Can you explain in a couple of paragraphs the different dynamic among these three countries?

The proportion of each part of the three countries are similar but with small differences. The gap between gross capital formation and household expenditure is larger when it comes to India, which means India government spends less money and the consumption of India is getting smaller and smaller. Gross capital formation and household expenditure are close in Germany an US.

If you want to, please change country_list <- c("United States","India", "Germany") to include your own country and compare it with any two other countries you like

We have chosen data of China, India and America. Comparing these three countries, the proportion of household expenditure in China is lower than the other countries and the proportion of gross capital formation is higher.

6 Deliverables

There is a lot of explanatory text, comments, etc. You do not need these, so delete them and produce a stand-alone document that you could share with someone. Knit the edited and completed R Markdown file as an HTML document (use the “Knit” button at the top of the script editor window) and upload it to Canvas.

7 Details

• Who did you collaborate with: Within our own group
• Approximately how much time did you spend on this problem set: Three days
• What, if anything, gave you the most trouble:

Please seek out help when you need it, and remember the 15-minute rule. You know enough R (and have enough examples of code from class and your readings) to be able to do this. If you get stuck, ask for help from others, post a question on Slack– and remember that I am here to help too!

As a true test to yourself, do you understand the code you submitted and are you able to explain it to someone else?

8 Rubric

Check minus (1/5): Displays minimal effort. Doesn’t complete all components. Code is poorly written and not documented. Uses the same type of plot for each graph, or doesn’t use plots appropriate for the variables being analyzed.

Check (3/5): Solid effort. Hits all the elements. No clear mistakes. Easy to follow (both the code and the output).

Check plus (5/5): Finished all components of the assignment correctly and addressed both challenges. Code is well-documented (both self-documented and with additional comments as necessary). Used tidyverse, instead of base R. Graphs and tables are properly labelled. Analysis is clear and easy to follow, either because graphs are labeled clearly or you’ve written additional text to describe how you interpret the output.